

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
lambdaStar = theta((size(thetaStar,1)-5):end,use);
UStar =  payoffs(thetaStar,X,Z,YY,reb,mpow);
clear K l0 srngSV SV theta theta0 time use Y exitn;


%% Preallocate
B = 2000;
PRpeace = NaN(B, N, 13);
Enumenter = NaN(B, N, 13);
Enumentercond = NaN(B, N, 13);
WW = NaN(B, N, 13);


%% expectations in data
load('sim_eq_final.mat')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);

parfor n=1:N    
    for b=1:B
    if ~isempty(intervEq{b,n})
        PrA = profprob(1:size(YY,1), DMES{b,n}, intervEq{b,n}, V_R(1,b,n), lambdaStar);
        PRpeace(b,n,1) = PrA(1,1);
        Enumenter(b,n,1) = PrA' * sum(YY(:,2:end),2);
        
        peq = exp(DMES{b,n} * lambdaStar);
        peq = peq / sum(peq);
        Enumentercond(b,n,1) = sum(sum((peq' .* intervEq{b,n}(2:end,:)) .* sum(YY(2:end,2:end),2)));
     end
    end
end


%% expectations in all/none intervene counterfactuals
Ustayout = squeeze(V_R(YY(:,1)==0,:,:));
Uenter_noInt = squeeze(V_R(YY(:,1)==1 & sum(YY(:,2:end),2)==0,:,:));
Uenter_allInt = squeeze(V_R(YY(:,1)==1 & sum(YY(:,2:end),2)==5,:,:));

PRpeace(:,:,2) = Ustayout >= Uenter_noInt;
PRpeace(:,:,3) = Ustayout >= Uenter_allInt;

Enumenter(:,:,2) =PRpeace(:,:,2)*5;
Enumenter(:,:,3) =PRpeace(:,:,3)*5;

clear Ustayout Uenter_noInt Uenter_allInt V_R V_UK V_France V_Russia V_China ans V_US DMES intervEq


%% expectations in US counterfactuals 
load('counterfactualUS_replication.mat')
ne4 = sum(YY4(:,2:end),2) + [0;ones(2^4,1)];

parfor n=1:N
    for b=1:B
        if ~isempty(noUSintervEq{b,n})
            PrAno =  profprob(1:size(YY4,1), noUSDMES{b,n}, noUSintervEq{b,n}, VnoUS_R(1,b,n), lambdaStar);
            PRpeace(b,n,4) =PrAno(1,1);
            Enumenter(b,n,4) = PrAno' * sum(YY4(:,2:end),2);
            
            peq = exp(noUSDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);
            Enumentercond(b,n,4) = sum(sum((peq' .* noUSintervEq{b,n}(2:end,:)) .* sum(YY4(2:end,2:end),2)));
        end
    end
end

parfor n=1:N
    for b=1:B
        if ~isempty(allUSintervEq{b,n})
            PrAall = profprob(1:size(YY4,1), allUSDMES{b,n}, allUSintervEq{b,n}, VallUS_R(1,b,n), lambdaStar);
            PRpeace(b,n,5) =PrAall(1,1);
            Enumenter(b,n,5) = PrAall' * ne4;
            
            peq = exp(allUSDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);
            Enumentercond(b,n,5) = sum(sum((peq' .* allUSintervEq{b,n}(2:end,:)) .* ne4(2:end)));
        end
    end
end

clear allUSDMES allUSeqPayoffs allUSintervEq ans noUSDMES noUSeqPayoffs noUSintervEq srng_cfactUS timeUS 
clear VallUS_China VallUS_France VallUS_R VallUS_Russia VallUS_UK VallUS_US 
clear VnoUS_China VnoUS_France VnoUS_R VnoUS_Russia VnoUS_UK VnoUS_US 



%% expectations in UK counterfactuals 
load('counterfactualUK_replication.mat')

parfor n=1:N
    for b=1:B
        if ~isempty(noUKintervEq{b,n})
            PrAno =  profprob(1:size(YY4,1), noUKDMES{b,n}, noUKintervEq{b,n}, VnoUK_R(1,b,n), lambdaStar);
            PRpeace(b,n,6) =PrAno(1,1);
            Enumenter(b,n,6) = PrAno' * sum(YY4(:,2:end),2);
            
            peq = exp(noUKDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);
            Enumentercond(b,n,6) = sum(sum((peq' .* noUKintervEq{b,n}(2:end,:)) .* sum(YY4(2:end,2:end),2)));
        end
    end
end

parfor n=1:N
    for b=1:B
        if ~isempty(allUKintervEq{b,n})
            PrAall = profprob(1:size(YY4,1), allUKDMES{b,n}, allUKintervEq{b,n}, VallUK_R(1,b,n), lambdaStar);
            PRpeace(b,n,7) =PrAall(1,1);
            Enumenter(b,n,7) = PrAall' * ne4;
            
            peq = exp(allUKDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);            
            Enumentercond(b,n,7) = sum(sum((peq' .* allUKintervEq{b,n}(2:end,:)) .* ne4(2:end)));
        end
    end
end
clear allUKDMES allUKeqPayoffs allUKintervEq ans noUKDMES noUKeqPayoffs noUKintervEq srng_cfactUK timeUK 
clear VallUK_China VallUK_France VallUK_R VallUK_Russia VallUK_UK VallUK_US 
clear VnoUK_China VnoUK_France VnoUK_R VnoUK_Russia VnoUK_UK VnoUK_US 



%% expectations in France counterfactuals 
load('counterfactualFR_replication.mat')

parfor n=1:N
    for b=1:B
        if ~isempty(noFRintervEq{b,n})
            PrAno =  profprob(1:size(YY4,1), noFRDMES{b,n}, noFRintervEq{b,n}, VnoFR_R(1,b,n), lambdaStar);
            PRpeace(b,n,8) =PrAno(1,1);
            Enumenter(b,n,8) = PrAno' * sum(YY4(:,2:end),2);
            
            peq = exp(noFRDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);            
            Enumentercond(b,n,8) = sum(sum((peq' .* noFRintervEq{b,n}(2:end,:)) .* sum(YY4(2:end,2:end),2)));
        end
    end
end

parfor n=1:N
    for b=1:B
        if ~isempty(allFRintervEq{b,n})
            PrAall = profprob(1:size(YY4,1), allFRDMES{b,n}, allFRintervEq{b,n}, VallFR_R(1,b,n), lambdaStar);
            PRpeace(b,n,9) =PrAall(1,1);
            Enumenter(b,n,9) = PrAall' * ne4;
            
            peq = exp(allFRDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);            
            Enumentercond(b,n,9) = sum(sum((peq' .* allFRintervEq{b,n}(2:end,:)) .* ne4(2:end)));        
        end
    end
end

clear allFRDMES allFReqPayoffs allFRintervEq ans noFRDMES noFReqPayoffs noFRintervEq srng_cfactFR timeFR 
clear VallFR_China VallFR_France VallFR_R VallFR_Russia VallFR_UK VallFR_US 
clear VnoFR_China VnoFR_France VnoFR_R VnoFR_Russia VnoFR_UK VnoFR_US 



%% expectations in Russia counterfactuals 
load('counterfactualRU_replication.mat')

parfor n=1:N
    for b=1:B
        if ~isempty(noRUintervEq{b,n})
            PrAno =  profprob(1:size(YY4,1), noRUDMES{b,n}, noRUintervEq{b,n}, VnoRU_R(1,b,n), lambdaStar);
            PRpeace(b,n,10) =PrAno(1,1);
            Enumenter(b,n,10) = PrAno' * sum(YY4(:,2:end),2);

            peq = exp(noRUDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);            
            Enumentercond(b,n,10) = sum(sum((peq' .* noRUintervEq{b,n}(2:end,:)) .* sum(YY4(2:end,2:end),2)));        
        end
    end
end

parfor n=1:N
    for b=1:B
        if ~isempty(allRUintervEq{b,n})
            PrAall = profprob(1:size(YY4,1), allRUDMES{b,n}, allRUintervEq{b,n}, VallRU_R(1,b,n), lambdaStar);
            PRpeace(b,n,11) =PrAall(1,1);
            Enumenter(b,n,11) = PrAall' * ne4;
            
            peq = exp(allRUDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);            
            Enumentercond(b,n,11) = sum(sum((peq' .* allRUintervEq{b,n}(2:end,:)) .* ne4(2:end)));  
        end
    end
end

clear allRUDMES allRUeqPayoffs allRUintervEq ans noRUDMES noRUeqPayoffs noRUintervEq srng_cfactRU timeRU 
clear VallRU_China VallRU_France VallRU_R VallRU_Russia VallRU_UK VallRU_US 
clear VnoRU_China VnoRU_France VnoRU_R VnoRU_Russia VnoRU_UK VnoRU_US



%% expectations in China counterfactuals 
load('counterfactualCH_replication.mat')

parfor n=1:N
    for b=1:B
        if ~isempty(noCHintervEq{b,n})
            PrAno =  profprob(1:size(YY4,1), noCHDMES{b,n}, noCHintervEq{b,n}, VnoCH_R(1,b,n), lambdaStar);
            PRpeace(b,n,12) =PrAno(1,1);
            Enumenter(b,n,12) = PrAno' * sum(YY4(:,2:end),2);
            
            peq = exp(noCHDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);            
            Enumentercond(b,n,12) = sum(sum((peq' .* noCHintervEq{b,n}(2:end,:)) .* sum(YY4(2:end,2:end),2)));                
        end
    end
end

parfor n=1:N
    for b=1:B
        if ~isempty(allCHintervEq{b,n})
            PrAall = profprob(1:size(YY4,1), allCHDMES{b,n}, allCHintervEq{b,n}, VallCH_R(1,b,n), lambdaStar);
            PRpeace(b,n,13) =PrAall(1,1);
            Enumenter(b,n,13) = PrAall' * ne4;
            
            peq = exp(allCHDMES{b,n} * lambdaStar);
            peq = peq / sum(peq);            
            Enumentercond(b,n,13) = sum(sum((peq' .* allCHintervEq{b,n}(2:end,:)) .* ne4(2:end)));
        end
    end
end
clear allCHDMES allCHeqPayoffs allCHintervEq ans noCHDMES noCHeqPayoffs noCHintervEq srng_cfactCH timeCH 
clear VallCH_China VallCH_France VallCH_R VallCH_Russia VallCH_UK VallCH_US 
clear VnoCH_China VnoCH_France VnoCH_R VnoCH_Russia VnoCH_UK VnoCH_US 



%% save to .csv

load('data.mat')
data.Properties.VariableNames{1,1} = 'ccode';
Enumentercond(isinf(Enumentercond))  = NaN;
meanPRpeace=NaN(N,size(PRpeace,3));
meanEnumenter = NaN(N,size(Enumenter,3));
meanEnumentercond = NaN(N,size(Enumentercond,3));

for i = 1:13
   meanPRpeace(:,i) = nanmean(PRpeace(:,:,i));
   meanEnumenter(:,i) = nanmean(Enumenter(:,:,i));
   meanEnumentercond(:,i) = nanmean(Enumentercond(:,:,i));
end

cf_names = {'data' 'noInt' 'allInt' 'noUS', 'allUS', ...
                'noUK', 'allUK', 'noFRN', 'allFRN', 'noRUS', 'allRUS', 'noCHN', 'allCHN'};
meanPRpeace = array2table(meanPRpeace, 'VariableNames', cf_names);
meanPRpeace = [data(:,1:2), meanPRpeace];

meanEnumenter = array2table(meanEnumenter, 'VariableNames', cf_names);
meanEnumenter = [data(:,1:2), meanEnumenter];

meanEnumentercond = array2table(meanEnumentercond, 'VariableNames', cf_names);
meanEnumentercond = [data(:,1:2), meanEnumentercond];

writetable(meanPRpeace, 'counterfactuals_prpeace_replication.csv')



